Analysis of First LIGO Science Data for Stochastic Gravitational Waves 
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We present the analysis of between 50 and 100 hrs of coincident interferometric strain data used 
to search for and estabUsh an upper limit on a stochastic background of gravitational radiation. 
These data come from the first LIGO science run, during which all three LIGO interferometers 
were operated over a 2-week period spanning August and September of 2002. The method of cross- 
correlating the outputs of two interferometers is used for analysis. We describe in detail practical 
signal processing issues that arise when working with real data, and we establish an observational 
upper limit on a f^^ power spectrum of gravitational waves. Our 90% confidence limit is Q,q /ifoo ^ 
23 in the frequency band 40 to 314 Hz, where ftioo is the Hubble constant in units of 100 km/sec/Mpc 
and fio is the gravitational wave energy density per logarithmic frequency interval in units of the 
closure density. This limit is approximately 10* times better than the previous, broadband direct 
limit using interferometric detectors, and nearly 3 times better than the best narrow-band bar 
detector limit. As LIGO and other worldwide detectors improve in sensitivity and attain their design 
goals, the analysis procedures described here should lead to stochastic background sensitivity levels 
of astrophysical interest. 
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I. INTRODUCTION 

In the last few years a number of new gravitational 
wave detectors, using long-baseline laser interferometry, 
have begun operation. These include the Laser Inter- 
ferometer Gravitational Wave Observatory (LIGO) de- 
tectors located in Hanford, WA and Livingston, LA Q; 
the GEO-600 detector near Hannover, Germany the 
VIRGO detector near Pisa, Italy Q; and the Japanese 
TAMA-300 detector in Tokyo While aU of these 

instruments are still being commissioned to perform at 
their designed sensitivity levels, many have begun mak- 
ing dedicated data collecting runs and performing gravi- 
tational wave search analyses on these data. 

In particular, from 23 August 2002 to 9 September 
2002, the LIGO Hanford and LIGO Livingston Observa- 
tories (LHO and LLO) collected coincident science data; 
this first scientific data run is referred to as SI. The 
LHO site contains two identically oriented interferome- 
ters: one having 4 km long measurement arms (referred 
to as HI), and one having 2 km long arms (H2); the 
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LLO site contains a single, 4 km long interferometer 
(LI). GEO-600 also took data in coincidence with the 
LIGO detectors during that time. Members of the LIGO 
Scientific Collaboration have been analyzing these data 
to search for gravitational wave signals. These initial 
analyses are aimed at developing the search techniques 
and machinery, and at using these fundamentally new in- 
struments to tighten upper limits on gravitational wave 
sources. Here we report on the methods and results of 
an analysis performed on the LIGO data to set an upper 
limit on a stochastic background of gravitational waves. ^ 
This represents the first such analysis performed on data 
from these new long-baseline detectors. The outline of 
the paper is as follows: 

Section II gives a description of the LIGO instruments 
and a summary of their operational characteristics during 
the SI data run. In Sec. Ill, we give a brief description of 
the properties of a stochastic background of gravitational 
radiation, and Sec. IV reviews the basic analysis method 
of cross-correlating the outputs of two gravitational wave 
detectors. 

In Sec. V, we discuss in detail the analysis performed 
on the LIGO data set. In applying the basic cross- 
correlation technique to real detector data, we have ad- 
dressed some practical issues and performed some ad- 
ditional analyses that have not been dealt with previ- 
ously in the literature: (i) avoidance of spectral leak- 
age in the short-time Fourier transforms of the data; (ii) 
a procedure for identifying and removing narrow-band 
(discrete frequency) correlations between detectors; (iii) 
chi-squared and time shift analyses, designed to explore 
the frequency domain character of the cross correlations. 

In Sec. VI, the error estimation is presented, and in 
Sec. VII, we show how the procedure has been tested by 
analyzing data that contain an artificially injected, sim- 
ulated stochastic background signal. Section VIII dis- 
cusses in more detail the instrumental correlation that is 
observed between the two Hanford interferometers (HI 
and H2), and Sec. IX concludes the paper with a brief 
summary and topics for future work. 

Appendix IXI gives a list of symbols used in the paper, 
along with their descriptions and equation numbers or 
sections in which they were defined. 



II. THE LIGO DETECTORS 

An interferometric gravitational-wave detector at- 
tempts to measure oscillations in the space-time metric, 
utilizing the apparent change in light travel time induced 
by a gravitational wave. At the core of each LIGO de- 
tector is an orthogonal arm Michelson laser interferom- 



Given the GEO SI sensitivity level and large geographical sepa- 
ration of the GEO-600 and LIGO detectors, it was not profitable 
to include GEO-600 data in this analysis. 



eter, as its geometry is well-matched to the space-time 
distortion. During any half-cycle of the oscillation, the 
quadrupolar gravitational-wave field increases the light 
travel time in one arm and decreases it in the other arm. 
Since the gravitational wave produces the equivalent of a 
strain in space, the travel time change is proportional to 
the arm length, hence the long arms. Each arm contains 
two test masses, a partially transmitting mirror near the 
beam splitter and a near-perfect reflector at the end of 
the arm. Each such pair is oriented to form a resonant 
Fabry-Perot cavity, which further increases the strain in- 
duced phase shifts by a factor proportional to the cav- 
ity finesse. An additional partially transmitting mirror 
is placed in the input path to form the power-recycling 
cavity, which increases the power incident on the beam 
splitter, thereby decreasing the shot-noise contribution to 
the signal-to-noisc ratio of the gravitational-wave signal. 

Each interferometer is illuminated with a medium 
power Nd:YAG laser, operating at 1.06 microns |^. Be- 
fore the light is launched into the interferometer, its fre- 
quency, amplitude and direction arc all stabilized, using 
a combination of active and passive stabilization tech- 
niques. To isolate the test masses and other optical ele- 
ments from ground and acoustic vibrations, the detectors 
implement a combination of passive and active seismic 
isolation systems 0, 0] , from which the mirrors are sus- 
pended as pendulums. This forms a coupled oscillator 
system with high isolation for frequencies above 40 Hz. 
The test masses, major optical components, vibration 
isolation systems, and main optical paths are all enclosed 
in a high vacuum system. 

Various feedback control systems arc used to keep the 
multiple optical cavities tightly on resonance Q and well- 
aligned [3- The gravitational wave strain signal is ob- 
tained from the error signal of the feedback loop used to 
control the differential motion of the interferometer arms. 
To calibrate the error signal, the effect of the feedback 
loop gain is measured and divided out, and the response 
R{ f) to a differential arm strain is measured and factored 
in. For the latter, the absolute scale is established us- 
ing the laser wavelength, and measuring the mirror drive 
signal required to move through a given number of in- 
terference fringes. During interferometer operation, the 
calibration was tracked by injecting fixed-amplitude sinu- 
soidal signals into the differential arm control loop, and 
monitoring the amplitude of these signals at the mea- 
surement (error) point jT3 |. 

Figure n shows reference amplitude spectra of equiva- 
lent strain noise, for the three LIGO interferometers dur- 
ing the SI run. The eventual strain noise goal is also in- 
dicated for comparison. The differences among the three 
spectra reflect differences in the operating parameters 
and hardware implementations of the three instruments; 
they are in various stages of reaching the flnal design 
configuration. For example, all interferometers operated 
during SI at a substantially lower effective laser power 
level than the eventual level of 6 W at the interferometer 
input; the resulting reduction in signal-to- noise ratio is 



4 

even greater than the square-root of the power reduction, 
because the detection scheme is designed to be efficient 
only near the design power level. Thus the shot-noise re- 
gion of the spectrum (above 200 Hz) is much higher than 
the design goal. Other major differences between the SI 
state and the final configuration were: partially imple- 
mented laser frequency and amplitude stabilization sys- 
tems; and partially implemented alignment control sys- 
tems. 




10^ Frequency (Hz) 10^ 



FIG. 1: Reference sensitivity curves for the three LIGO in- 
terferometers during the SI data run, in terms of equiva- 
lent strain noise density. The HI and H2 spectra are from 
9 September 2002, and the LI spectrum is from 7 Septem- 
ber 2002. Also shown are strain spectra corresponding to two 
levels of a stochastic background of gravitational radiation de- 
fined by Eq.EIl These can be compared to the expected 90% 
confidence level upper limits, assuming Gaussian uncorrelated 
detector noise at the levels shown here, for the interferome- 
ter pairs: H2-L1 (fio /iioo = 10), with 100 hours of correlated 
integration time; H1-H2 (fio /ifgo — 0.83), with 150 hours of 
integration time; and Hl-Ll (ilo ftioo = 11), with 100 hours 
of integration time. Also shown is the strain noise goal for 
the two 4-km arm interferometers (HI and LI). 

Two other important characteristics of the instru- 
ments' performance are the stationarity of the noise, and 
the duty cycle of operation. The noise was significantly 
non-stationary, due to the partial stabilization and con- 
trols mentioned above. In the frequency band of most 
importance to this analysis, approximately 60 — 300 Hz, 
a factor of 2 variation in the noise amplitude over several 
hours was typical for the instruments; this is addressed 
quantitatively in Sec. IVII and Fig. ^| As our analysis 
relies on cross-correlating the outputs of two detectors, 
the relevant duty cycle measures are those for double- 
coincident operation. For the SI run, the total times of 
coincident science data for the three pairs are: H1-H2, 
188 hrs (46% duty cycle over the SI duration); Hl-Ll, 
116 hrs (28%); H2-L1, 131 hrs (32%). A more detailed 
description of the LIGO interferometers and their per- 
formance during the SI run can be found in reference 
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III. 



STOCHASTIC GRAVITATIONAL WAVE 
BACKGROUNDS 

A. Spectrum 



A stochastic background of gravitational radiation is 
analogous to the cosmic microwave background radi- 
ation, though its spectrum is unlikely to be thermal. 
Sources of a stochastic background could be cosmologi- 
cal or astrophysical in origin. Examples of the former are 
zero-point fluctuations of the space-time metric amplified 
during inflation, and first-order phase transitions and de- 
caying cosmic string networks in the early universe. An 
example of an astrophysical source is the random super- 
position of many weak signals from binary-star systems. 
See references [TjI and [13 for a review of sources. 

The spectrum of a stochastic background is usually 
described by the dimensionless quantity f2gw(/) which is 
the gravitational-wave energy density per unit logarith- 
mic frequency, divided by the critical energy density pc 
to close the universe: 



Pc df 



(3.1) 



The critical density pc = 3c Hq /8nG depends on the 
present day Hubble expansion rate Hq. For convenience 
we define a dimensionless factor 



where 

-ffioo = 100 



km 



Ho/H_ 



100 : 



(3.2) 



sec • Mpc 



3.24x10-1^ — , (3.3) 
sec 



to account for the different values of Hq that are quoted 
in the literature.^ Note that rigw(/) /iioo is independent 
of the actual Hubble expansion rate, so we work with this 
quantity rather than Qg^{f) alone. 

Our specific interest is the measurable one-sided power 
spectrum of the gravitational wave strain S'g„(/), which 
is normalized according to: 



hm 1 / ^ dt\hit)\' 



dfSg^if), (3.4) 



where h(t) is the strain in a single detector due to the 
gravitational wave signal; h{t) can be expressed in terms 



2 ifo = 73±2±7 km/sec/Mpc as shown in Ref. Q and from inde- 
pendent SNIa observations from observatories on the ground Il5l . 
The Wilkinson Microwave Anisotropy Pro be 1st year (WMAPl) 
observation has Hq = Tl^g km/sec/Mpc Il6l . 



of the perturbations hab of the spacetime metric and the 
detector geometry via: 



h{t)= habit, Xo)^{X''X' 



(3.5) 



Here xq specifies the coordinates of the interferometer 
vertex, and y° are unit vectors pointing in the di- 
rection of the detector arms. Since the energy density in 
gravitational waves involves a product of time derivatives 
of the metric perturbations (c.f. p. 955 of Ref. Il7['). one 
can show (see, e.g.. Sees. H.A and HI. A in Ref. [I3 for 
more details) that Sgw{f) is related to ilgw(/) via: 



5'gw(/) — 



^■^0 f-3n (f\ 



(3.6) 



Thus, for a stochastic gravitational wave background 
with rigw(/) = f^o = const (as is predicted at LIGO 
frequencies e.g., by inflationary models in the infinitely 
slow- roll limit; or by cosmic string models |l9l |') the power 
in gravitational waves falls off as 1//^, with a strain am- 
plitude scale of: 



5.6 X 10 



-22 



noo- 



i2^)'^'Hz-i/2. (3.7) 



The spectrum f2gw(/) completely specifies the statis- 
tical properties of a stochastic background of gravita- 
tional radiation provided we make several additional as- 
sumptions. Here, we assume that the stochastic back- 
ground is isotropic, unpolarized, stationary, and Gaus- 
sian. Anisotropic or non-Gaussian backgrounds (e.g., 
due to an incoherent superposition of gravitational waves 
from a large number of unresolved white dwarf binary 
star systems in our own galaxy, or a "pop-corn" stochas- 
tic signal produced by gravitational waves from super- 
nova core-collapse events pol \2l\) may require different 
data analysis techniques from those presented here. (See, 
e-g-, [2^ 113 for a discussions of these different tech- 
niques.) 



B. Prior observational constraints 

While predictions for rigw(/) from cosmological mod- 
els can vary over many orders of magnitude, there are 
several observational results that place interesting upper 
limits on rig„(/) in various frequency bands. Table HI 
summarizes these observational constraints and upper 
limits on the energy density of a stochastic gravitational 
wave background. The high degree of isotropy observed 
in the cosmic microwave background radiation (CMBR) 
places a strong constraint on r2gw(/) at very low frequen- 
cies |23|- Since Hiqq sa 3.24 x 10"^^ Hz, this hmit apphes 
only over several decades of frequency 10~^* — 10~^^ Hz 
which are far below the bands accessible to investiga- 
tion by either Earth-based (10 — 10"* Hz) or space-based 
(10-* - 10-1 Hz) detectors. 
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Another observational constraint comes from nearly 
two decades of monitoring the time-of-arrival jitter of 
radio pulses from a number of millisecond pulsars [25l |. 
These pulsars are remarkably stable clocks, and the reg- 
ularity of their pulses places tight constraints on Ogw(/) 
at frequencies on the order of the inverse of the obser- 
vation time of the pulsars, 1/T ~ 10~^ Hz. Like the 
constraint derived from the isotropy of the CMBR, the 
millisecond pulsar timing constraint applies to an obser- 
vational frequency band much lower than that probed by 
Earth-based and space-based detectors. 

The only constraint on f2gw(/) within the frequency 
band of Earth-based detectors comes from the observed 
abundances of the light elements in the universe, cou- 
pled with the standard model of big-bang nucleosynthe- 
sis . For a narrow range of key cosmological param- 
eters, this model is in remarkable agreement with the 
elemental observations. One of the constrained parame- 
ters is the expansion rate of the universe at the time of 
nucleosynthesis, thus constraining the energy density of 
the universe at that time. This in turn constrains the 
energy density in a cosmological background of gravita- 
tional radiation (non-cosmological sources of a stochastic 
background, e.g., from a superposition of supernovae sig- 
nals, are not of course constrained by these observations). 
The observational constraint is on the logarithmic inte- 
gral over frequency of f7gw(/). 

All the above constraints were indirectly inferred via 
electromagnetic observations. There are a few, much 
weaker constraints on ^gw{f) that have been set by ob- 
servations with detectors directly sensitive to gravita- 
tional waves. The earliest such measurement was made 
with room-temperature bar detectors, using a split bar 
technique for wide bandwidth performance p^ . Later 
measurements include an upper limit from a correlation 
between the Garching and Glasgow prototype interfer- 
ometers , several upper limits from observations with 
a single cryogenic resonant bar detector [2^ l30l |. and 
most recently an upper limit from observations of two- 
detector correlations between the Ex plor er and Nautilus 
cryogenic resonant bar detectors |3l|, |32| ■ Note that the 
cryogenic resonant bar observations are constrained to a 
very narrow bandwidth (Af ~ IHz) around the resonant 
frequency of the bar. 

IV. DETECTION VIA CROSS-CORRELATION 

We can express the equivalent strain output Si{t) of 
each of our detectors as: 

s,{t) = h,(t) + n,{t) , (4.1) 

where hi{t) is the strain signal in the i-th detector due 
to a gravitational wave background, and Ui (i) is the de- 
tector's equivalent strain noise. If we had only one de- 
tector, all we could do would be to put an upper limit 
on a stochastic background at the detector's strain noise 
level; e.g., using LI wc could put a limit of fig /i^oo ~ 10^ 



in the band 100 — 200 Hz. To do much better, we cross- 
correlate the outputs of two detectors, taking advantage 
of the fact that the sources of noise Uiin each detector 
will, in general, be independent [H El [H IH H HI ■ 
We thus compute the general cross-correlation:'^ 

Y = [ dh [ dt2Si{ti)Q{h-t2)s2{t2), (4.2) 

J-T/2 J-T/2 

where Q{ti — is a (real) filter function, which wc will 
choose to maximize the signal-to-noise ratio of Y. Since 
the optimal choice of Q{ti — ^2) fahs off rapidly for time 
delays \ti — ^2! large compared to the light travel time 
d/c between the two detectors,'' and since a typical ob- 
servation time T will be much, much greater than d/c, 
we can change the limits on one of the integrations from 
(— T/2,r/2) to (—00,00), and subsequently obtain [r§ |: 

df / df dT{f-f')~si{f)Q{r)s2ir), 

-00 'J —00 

(4.3) 

where 

Srif)^ r\te---f^='^^ (4.4) 

is a finite-time approximation to the Dirac delta function, 
and Si(/), Q{f) denote the Fourier transforms of Si{t), 
Q(t)-i.e., aif) EE e-2./t 

To find the optimal Q{f), we assume that the intrin- 
sic detector noise is: (i) stationary over a measurement 
time T; (ii) Gaussian; (iii) uncorrelated between different 
detectors; (iv) uncorrelated with the stochastic gravita- 
tional wave signal; and (v) much greater in power at any 
frequency than the stochastic gravitational wave back- 
ground. Then the expected value of the cross-correlation 
Y depends only on the stochastic signal: 

fiY^{Y)^-J d/7(l/l)^gw(|/|)Q(/), (4.5) 

while the variance of Y is dominated by the noise in the 
individual detectors: 

4^((r-(y))2)«- / d/Fi(|/|)|Q(/)pP2(|/|). 

^ J — 00 

(4.6) 

Here Pi{f) and P2{f) are the one-sided strain noise 
power spectra of the two detectors. The integrand of 
Eg. 14.51 contains a (real) function 7(/), called the over- 
lap reduction function |35l | . which characterizes the re- 
duction in sensitivity to a stochastic background arising 



^ The equ ations in this section are a summary of Sec. Ill from 
Ref. 11^ . Readers interested in mo re d etails and/or derivations 
of the key equations should refer to llSi and references contained 
therein. 

* The light travel time d/c between the Hanford and Livingston 
detectors is approximately 10 msec. 
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Observational 
Technique 


Observed 
Limit 


Frequency 
Domain 


Comments 


Cosmic Microwave 
Background 


f^gw(/) fe?oo<10-^-^(ifi^)' 


3 X 10~^**Hz < / < 10"^'' Hz 


[24] 


Radio Pulsar 
Timing 


^^gw(/) hloo < 9.3 X 10"* 


4 X 10"^ Hz < / < 4 X 10"* Hz 


95% CL bound, [25j 


Big-Bang 
Nucleosynthesis 


//>io-B HzC'ln/ ^^gw(/) hloo < 10-' 


/ > 10"*Hz 


95% CL bound, [2^ 


Interferometers 


f^gw(/) hloo < 3 X 10^ 


100 Hz < f < 1000 Hz 


Garching-Glasgow 


Room Temperature 
Resonant Bar (correlation) 


!^gw(/o) hloo < 3000 


fo = 985 ± 80 Hz 


Clasgow [27] 


Cryogenic Resonant 
Bar (single) 


f^gw(/o) ft?oo <300 
f^gw(/o) hloo < 5000 


fo = 907 Hz 
fo = 1875 Hz 


Explorer [29j 
ALTAIR [30] 


Cryogenic Resonant 
Bar (correlation) 


!^gw(/o) hloo < 60 


fo = 907 Hz 


Explorer+Nautilus [31. 321 



TABLE I: Summary of upper limits on hloo o^^i' ^ large range of frequency bands. The upper portion of the table lists 
indirect limits derived from astrophysical observations. The lower portion of the table lists limits obtained from prior direct 
gravitational wave measurement. 



from the separation time delay and relative orientation of 
the two detectors. It is a function of only the relative de- 
tector geometry (for coincident and co- aligned detectors, 
like HI and H2, j{f) = 1 for all frequencies). A plot of 
the overlap reduction function for correlations between 
LIGO Livingston and LIGO Hanford is shown in Fig. |2 



From Eqs. 14.511 and (|4.t)|l . it is relatively straightfor- 
ward to show |l2| that the expected signal-to-noise ratio 
(/iy/(Ty) of Y is maximized when 



Qif) 



7(l/l)ggw(|/|) 

^(1/1)^2(1/1) 



(4.7) 



C--0.4 




150 200 250 

Frequency (Hz) 



For the SI analysis, we specialize to the case figw(f) 
flo = const. Then, 



Qif) =AA 



7(1/1) 



|/PPl(|/|)P2(|/|) 



(4.8) 



where A/" is a (real) overall normalization constant. 
In practice we choose J\f so that the expected cross- 



correlation is /iy = Qa /iJqo T. For such a choice. 



•iff 2 



df 



7^(1/1) 



.f'Pl{\f\)P2i\.f\) 



T 



IOtt^ 

3^100 



df 



7^(1/1) 



.m(i/)P2(i/i) 



(4.9) 

-1 

(4.10) 



In the sense that Q{f) maximizes ^.y/uy, it is the opti- 
mal filter for the cross-correlation Y . The signal-to-noise 
ratio py = Y jay has expected value 



FIG. 2: Overlap reduction function between the LIGO Liv- 
ingston and the LIGO Hanford sites. The value of I7I is a little 
less than unity at Hz because the interferometer arms are 
not exactly co-planar and co-aligned between the two sites. 



cry 

IOtt^ 



df 



iWf\) 



f^pm)mf\) 



(4.11) 

1/2 
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which grows with the square-root of the observation time 
T, and inversely with the product of the amplitude noise 
spectral densities of the two detectors. In order of magni- 
tude, Eg . 14 . 1 ll indicates that the upper limit we can place 
on riQ hiQQ by cross-correlation is smaller (i.e., more con- 
straining) than that obtainable from one detector by a 
factor of 7rms\/2^^BW5 where Abw is the bandwidth over 
which the integrand of Eg. IT. 1 H is significant (roughly the 
width of the peak of 1//^P,(/)), and Trms is the rms value 
of 7(/) over that bandwidth. For the LHO-LLO correla- 
tions in this analysis, T ^ 2 x 10^ sec, Abw ^ Hz, 
and 7rms ~ 0.1, so we expect to be able to set a limit that 
is a factor of several hundred below the individual detec- 
tors' strain noise^, or ilo ^loo ^ 10 as shown in Fig. ^ 



Detector #1 data 



Detector #2 data 



Dedmate to 
1024 samples/sec 



/ software \ 
injections ' 



Dedmate to 
1024 samples/sec 



Divide into 
1 segments 



Form power 
spedrum estimSe 
(0.25 Hz 



Form power 
spedrum estimde 
(0.25 Hz 



Divide into 
1 segments 



Apply calibration 



Apply calibBtion 



Conpute optimal 
filter & itieoretical 
variancefjv, ,2 



i. apply windov 

ii. ^ro pad 

iii. perform FFT 



Conpute ctoss- 
correlation Y| j 



i. apply windcw 

ii. ^ro pad 

iii. perform FFT 



These computations^formed on each 90-sec segme, )t 



Conpute interval 

mean Y\ & 
std da/iation s^^ 




FIG. 3: Data analysis flow diagram for the stochastic search. 
The raw detector signal (i.e., the uncalibrated differential arm 
error signal) is fed into the pipeline in 900-sec long intervals. 
Simulated stochastic background signals can be injected near 
the beginning of each data path, allowing us to test the data 
analysis routines in the presence of known correlations. 



ANALYSIS OF LIGO DATA 



A. Data analysis pipeline 



A flow diagram of the data analysis pipeline is shown 
in Fig. 131 \^(^ ■ We perform the analysis in the frequency 



^ More precisely, if the two detectors have unequal strain sensitiv- 
ities, the cross-correlation limit will be a factor of 7rms\/^'^BW 
below the geometric mean of the two noise spectral densities. 



domain, where it is more convenient to construct and 
apply the optimal filter. Since the data are discretely- 
sampled, we use discrete Fourier transforms and sums 
over frequency bins rather than integrals. The data ri[k] 
are the raw (uncalibrated) detector outputs at discrete 
times tk = k St: 



ri[k] = ri{tk) , 



(5.1) 



where k = 0, 1, 2, • • • , 5t is the sampling period, and i 
labels the detector. We decimate the data to a sam- 
pling rate of {5t)-'^ = 1024 Hz (from 16384 Hz), since 
the higher frequencies make a negligible contribution to 
the cross-correlation. The decimation is performed with 
a finite impulse response filter of length 320, and cut- 
off frequency 512 Hz. The data are split into intervals 
(labeled by index /) and segments (labeled by index J) 
within each interval to deal with detector nonstationar- 
ity and to produce sets of cross-correlation values Yjj for 
which empirical variances can be calculated; see Fig. 0] 
The time-series data corresponding to the J-th segment 
in interval / is denoted riij[k], where fc = 0, 1,-- - ,N — 1 
runs over the total number of samples in the segment. 

A single optimal filter Qi is calculated and applied 
for each interval /, the duration of which should be long 
enough to capture relatively narrow-band features in the 
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FIG. 4: Time-series data from each interferometer is split 
into Af 900-sec intervals, which are further subdivided into 
n = 10 90-sec data segments. Cross-correlation values Yjj are 
calculated for each 90-sec segment; theoretical variances o-y^j 
are calculated for each 900-sec interval. Here I — 1,2, ■ • ■ ,M 
labels the different intervals, and J — 1,2, •■• ,n labels the 
individual segments within each interval. 



power spectra, yet short enough to account for significant 
non-stationary detector noise. Based on observations of 
detector noise variation, we chose an interval duration of 
Tint = 900 sec. The segment duration sliould be much 
greater than the light travel time between the two de- 
tectors, yet short enough to yield a sufficient number of 
cross-correlation measurements within each interval to 
obtain an experimental estimate of the theoretical vari- 
ance J of the cross correlation statistic Yjj. We chose 
a segment duration of Tgcg = 90 sec, yielding ten Yjj 
values per interval. 

To compute the segment cross-correlation values Yjj, 
the raw decimated data rnj [k] are windowed in the time 
domain (see Sec. IV Bl for details), zero-padded to twice 
their length (to avoid wrap-around problems |37j when 
calculating the cross-correlation statistic in the frequency 
domain), and discrete Fourier-transformed. Explicitly, 
defining 



9iij[ 



^ I w,[k]nij[k] k^O,--- ,N-1 
I k^N,--- ,2N-1. 



(5.2) 



where Wi[k] is the window function for the i-th detector,^ 
the discrete Fourier transform is: 



2N-1 



g,ij[q] = J2 9uj[k]e-'^^ 



(5.3) 



frequency resolution, 5/, of the optimal filter Qf? 



1 



9*iiA<i]92ij[q] 

q—ribi—m 



(5.4) 



where ^min < ^ < ^max, = 2TscgSf is the number 
of frequency values being binned, and m = (rib — l)/2. 
The index £ labels the discrete frequencies, fe = £Sf. 
The are computed for a range of £ that includes 

only the frequency band that yields most of the expected 
signal-to-noise ratio (e.g., 40 to 314 Hz for the LHO- 
LLO correlations), as described in Sec. IV CI The cross- 
correlation values are calculated as: 



ij 



2 Re 



£ Sf Qi[£]gij[£] 



(5.5) 



Some of the frequency bins within the {£min, ^max} range 
are excluded from the above sum to avoid narrow-band 
instrumental correlations, as described in Sec. IV CI Ex- 
cept for the details of windowing, binning, and band- 
limiting, Eq. 15.51 for Yjj is just a discrete-frequency ap- 
proximation to Ea. l4.3l for Y for the continuous- frequency 
data, with df StU — /') approximated by a Kronecker 
delta 5ui in discrete frequencies and /^.^ 

In calculating the optimal filter, we estimate the strain 
noise power spectra Pu for the interval / using Welch's 
method: 449 periodograms are formed and averaged from 
4096-point, Hann- windowed data segments, overlapped 
by 50%, giving a frequency resolution Sf = 0.25 Hz. 
To calibrate the spectra in strain, we apply the cali- 
bration response function Ri{f) which converts the raw 
data to equivalent strain: Jiif) = R~^{f)ri{f). The cal- 
ibration lines described in Sec. ^] were measured once 
per 60 seconds; for each interval /, we apply the re- 
sponse function, Ru , corresponding to the middle 60 sec- 
onds of the interval. The optimal filter Qi for the case 
f^gw(/) = SIq = const is then constructed as: 



Qi 



\fi\ 



iRu[e]Pu[£])*{R: 



(5.6) 



2/ 



2/ 



fc=0 



where j[£] = jifi), and R^i[£] = Ruifi)-^By including 
the additional response function factors Ru in Eq. 15.61 
Qi has the appropriate units to act directly on the raw 
detector outputs in the calculation of Yjj (c.f. Ea. l5.5|) . 

The normalization factor A fr in Eq. 15.61 takes into ac- 
count the effect of windowing [33 . Choosing A// so that 



where N = Tgcg/St ~ 92160 is the number of data points 
in a segment, and q ~ 0, 1, • • • ,2N — 1. The cross- 
spectrum gljj[q] ■ ff2/,/['?] is formed and binned to the 



In general, one can use different window functions for different 
detectors. However, for the SI analysis, we took -ujiffe] = «)2[fc]. 



^ As discussed below, 5f = 0.25 Hz yielding nj, = 45 and m = 22. 

* To make this correspondence with Eq. 14.31 the factor of 2 and 
real part in Eq. 15.51 are needed since we are summing only over 
positive frequencies, e.g., 40 to 314 Hz for the LHO-LLO correla- 
tion. Basically, integrals over continuous frequency are replaced 
by sums over discrete frequency bins using the correspondence 
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the theoretical mean of the cross-correlation Yjj is equal 
to flo h^QQ Tscg for all I, J (as was done for Y in Sec. llV|l . 
we have: 



^fI 



1 



(JYij — Tscg 



20^2 



3i7? wiW2 



2 E ^/ 

IOtt 



(5.7) 



7 



2 X 2 



,,2, ,,2 



2 E ^/ 



(W1W2)^ 

2r/?i 



(5.8) 



7 



where 



N-l 



Wl'W2 



fc=0 

7V-1 



1^2 ^ ^E«^lW^2[fc] 



(5.9) 
(5.10) 



fe=0 



provided the windowing is sufRcient to prevent significant 
leakage of power across the frequency band (see Sec. IV Bl 
and for more details). Note that the theoretical vari- 
ance cTy^-^ depends only on the interval /, since the cross- 
correlations Yjj have the same statistical properties for 
each segment J in /. 

For each interval /. we calculate the mean, Yj, and 
(sample) standard deviation, syu, of the 10 cross- 
correlation values Yjj: 



1 ^" 
Yi = -TYjj. 
10^ 



J=i 



1 



(5.11) 
(5.12) 



We also form a weighted average, Kopt , of the Yj over the 
whole run: 



Yr 



opt 



(5.13) 



The statistic Y^pt maximizes the expected signal-to-noise 
ratio for a stochastic signal, allowing for non-statio nary 
detector noise from one 900-sec interval / to the next [l^ . 
Dividing Kopt by the time Tgog over which an individ- 
ual cross-correlation measurement is made gives, in the 
absence of cross-correlated detector noise, an unbiased 
estimate of the stochastic background level:^ /i^qo = 

Yopt/-^scg- 



^ Wc use a hat ^ to indicate an estimate of the actual (unknown) 
value of a quantity. 



Finally, in Sec. IV E( we will be interested in the spectral 
properties of Y/j, Yj, and i^pt. Thus, for later reference, 
we define: 



Yij[£] = QMGij 



Yj 



i^i:^"i 



./=i 



Y, 



opt [ 



(5.14) 
(5.15) 

(5.16) 



Note that 2 Re Ef^F i^/ ■ of the above quantities equal 
Y/j, Yf, and Yopt, respectively. 



B. Windowing 

In taking the discrete Fourier transform of the raw 90- 
sec data segments, care must be taken to limit the spec- 
tral leakage of large, low-frequency components into the 
sensitive band. In general, some combination of high- 
pass filtering in the time domain, and windowing prior 
to the Fourier transform can be used to deal with spec- 
tral leakage. In this analysis we have found it sufficient 
to apply an appropriate window to the data. 

Examining the dynamic range of the data helps es- 
tablish the allowed leakage. Figure ^ shows that the 
lowest instrument noise around 60 Hz is approximately 
10^^^/\/Hz (for LI). While not shown in this plot, the 
rms level of the raw data corresponds to a strain of or- 
der 10^^^, and is due to fluctuations in the 10 — 30 Hz 
band. Leakage of these low-frequency components must 
be at least below the sensitive band noise level; e.g., leak- 
age must be below lO^'' for a 30 Hz offset. A tighter 
constraint on the leakage comes when considering that 
these low-frequency components may be correlated be- 
tween the two detectors, as they surely will be at some 
frequencies for the two interferometers at LHO, due to 
the common seismic environment. In this case the leak- 
age should be below the predicted stochastic background 
sensitivity level, which is approximately 2.5 orders of 
magnitude below the individual detector noise levels for 
the LHO H1-H2 case. Thus, the leakage should be below 
3 X 10"^ for a 30 Hz offset. 

On the other hand, we prefer not to use a window 
that has an average value significantly less than unity 
(and correspondingly low leakage, such as a Hann win- 
dow), because it will effectively reduce the amount of 
data contributing to the cross-correlation. Provided that 
the windowing is sufficient to prevent significant leakage 
of power across the frequency range, the net effect is to 
multiply the expected value of the signal-to-noise ratio 

by wiW2/\Jwlwl (c.f. Eas. lHnO]) . 

For example, when wi and 1x12 are both Hann windows, 
this factor is equal to a/18/35 ~ 0.717, which is equiva- 
lent to reducing the data set length by a factor of 2. In 
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principle one should be able to use overlapping data seg- 
ments to avoid this effective loss of data, as in Welch's 
power spectrum estimation method. In this case, the 
calculations for the expected mean and variance of the 
cross-correlations would have to take into account the 
statistical interdependence of the overlappin g da ta. 

Instead, we have used a Tukcy window |39l| . which 
is essentially a Hann window split in half, with a con- 
stant section of all I's in the middle. We can choose the 
length of the Hann portion of the window to provide suffi- 
ciently low leakage, yet maintain a unity value over most 
of the window. FigureElshows the leakage function of the 
Tukey window that we use (a 1-sec Hann window with an 
89-sec flat section spliced into the middle) , and compares 
it to Hann and rectangular windows. The Tukey window 
leakage is less than 10~^ for all frequencies greater than 
35 Hz away from the FFT bin center. This is 4 orders of 
magnitude better than what is needed for the LHO-LLO 
correlations and a factor of 30 better suppression than 
needed for the H1-H2 correlation. 




Offset from FFT bin center (Hz) 



FIG. 5: Leakage function for a rectangular window, a stan- 
dard Hann window of width 90-sec, and a Tukey window con- 
sisting of a 1-sec Hann window with an 89-sec flat section 
spliced into the middle. The curves show the envelope of the 
leakage functions, with a varying frequency resolution, so the 
zeros of the functions are not seen. 

To explicitly verify that the Tukey window behaved 
as expected, we re-analyzed the H1-H2 data with a pure 
Hann window (see also Sec. IVIIIII . The result of this 
re-analysis, properly scaled to take into account the ef- 
fective reduction in observation time, was, within error, 
the same as the original analysis with a Tukey window. 
Since the H1-H2 correlation is the most prone of all corre- 
lations to spectral leakage (due to the likelihood of cross- 
correlated low- frequency noise components) , the lack of a 
significant difference between the pure Hann and Tukey 
window analyses provided additional support for the use 
of the Tukcy window. 



C. Frequency band selection and discrete 
frequency elimination 

In computing the discrete cross-correlation integral, we 
are free to restrict the sum to a chosen frequency region 
or regions; in this way the variance can be reduced (e.g., 
by excluding low frequencies where the detector power 
spectra arc large and relatively less stationary), while 
still retaining most of the signal. We choose the fre- 
quency ranges by determining the band that contributes 
most of the expected signal-to-noise ratio, according to 
Eq. 14.111 Using the strain power spectra shown in Fig. ^ 
we compute the signal-to-noise ratio integral of Eq. 14.111 
from a very low frequency (a few Hz) up to a variable 
cut-off frequency, and plot the resulting signal-to-noise 
ratio versus cut-off frequency (Fig. 1^. For each interfer- 
ometer pair, the lower band edge is chosen to be 40 Hz, 
while the upper band edge choices are 314 Hz for LHO- 
LLO correlations (where there is a zero in the overlap 
reduction function), and 300 Hz for H1-H2 correlations 
(chosen to exclude ~ 340 Hz resonances in the test mass 
mechanical suspensions, which were not well-resolved in 
the power spectra). 
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FIG. 6; Curves show the fraction of maximum expected 
signal-to-noise ratio as a function of cut-off frequency, for the 
three interferometer pairs. The curves were made by numer- 
ically integrating Eq. 14.111 from a few Hz up to the variable 
cut-off frequency, using the strain sensitivity spectra shown 
in Fig. [2 

Within the 40 - 314 (300) Hz band, discrete fre- 
quency bins at which there are known or potential in- 
strumental correlations due to common periodic sources 
arc eliminated from the cross-correlation sum. For ex- 
ample, a significant feature in all interferometer out- 
puts is a set of spectral lines extending out to beyond 
2 kHz, corresponding to the 60 Hz power line and its 
harmonics (n • 60 Hz lines). Since these lines obvi- 
ously have a common source — the mains power supplying 
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the instrumentation — they are potentially correlated be- 
tween detectors. To avoid including any such correlation 
in the analysis, we eliminate the n • 60 Hz frequency bins 
from the sum in Eq. 15.51 

Another common periodic signal arises from the data 
acquisition timing systems in the detectors. The abso- 
lute timing and synchronization of the data acquisition 
systems between detectors is based on 1 pulse-per-second 
signals produced by Global Positioning System (GPS) re- 
ceivers at each site. In each detector, data samples are 
stored temporarily in 1/16 sec buffers, prior to being col- 
lected and written to disk. The process, through mecha- 
nisms not yet established, results in some power at 16 Hz 
and harmonics in the detectors' output data channels. 
These signals are extremely narrow-band and, due to the 
stability and common source of the GPS-derivcd timing, 
can be correlated between detectors. To avoid including 
any of these narrow-band correlations, we eliminate the 
n ■ 16 Hz frequency bins from the sum in Eq. 15.51 

Finally, there may be additional correlated narrow- 
band features due to highly stable clocks or oscillators 
that are common components among the detectors (e.g., 
computer monitors can have very stable sync rates, typi- 
cally at 70 Hz). To describe how we avoid such features, 
we first present a quantitative analysis of the effect of 
coherent spectral lines on our cross-correlation measure- 
ment. We begin by following the treatment of correlated 
detector noise given in Sec. V.E of Ref. 0|. The con- 
tribution of cross-correlated detector noise to the cross- 
correlation Y will be small compared to the intrinsic mea- 
surement noise if 



< cry 



(5.17) 



where /'i2(|/|) is the cross-power spectrum of the strain 
noise (ni, 712) in the two detectors, T is the total obser- 
vation time, and ay is defined by Eg. 14.61 Using Eq. 
this condition becomes 



df 



^2(1/1)7(1/1) 



|/PPl(|/|)P2(|/|) 



2aY 



or, equivalently. 



57r2 



df 



^12(1/1)7(1/1) 



|/PPl(|/|)P2(|/|) 



< 2ct: 



(5.18) 



(5.19) 



where Eas. 14. 9| 14.101 were used to eliminate M in terms 



of ( 



57r2 



(5.20) 



Now consider the presence of a correlated periodic signal, 
such that the cross-spectrum Pi2{f) is significant only 
at a single (positive) discrete frequency, /l- For this 
component to have a small effect, the above condition 
becomes: 



57r2 



A/ 



Pi2(/l)7(/l) 



flPl{fL)P2{lL) 



(5.21) 



where A/ is the frequency resolution of the discrete 
Fourier transform used to approximate the frequency in- 
tegrals, and the left-hand-side of Eq. 15.211 can be ex- 
pressed in terms of the coherence function ri2(/), which 
is essentially a normalized cross-spectrum, defined as ^3 ■ 



^ \P12U)? 

^ 12{J) = 



The condition on the coherence at /l is thus 



w (f^^^n^^''y V^i(/l)P2(/l) 

[1 12(/l)J < -TJ |j_3 ,j, s| 



A/ iH^ 



100 



I/Z'7(/l)I 



(5.22) 



(5.23) 



Since cry increases as T^^"^, the limit on the coher- 
ence ri2(/L) becomes smaller as 1/T. To show how this 
condition applies to the SI data, we estimate the fac- 
tors in Ea. l5.23| for the H2-L1 pair, focusing on the band 
100 — 150 Hz. We assume any correlated spectral line 
is weak enough that it does not appear in the power 
spectrum estimates used to construct the optimal filter. 
Noting that the combination (STJ^^q/IOtt^)/"'^ is just 
the power spectrum of gravitational waves iS'gw(/) with 
^0 ^100 = 1 (*^-f- Ea- I3.6|l . we can evaluate the right-hand- 
side of Eq.|^23by estimating the ratios [f'i/S'gw]^^^ from 
Fig. n for fio /i?oo = 1- Within the band 100 - 150 Hz, 
this gives: (PiP2)^^^/5'gw > 2500. The overlap reduction 
function in this band is I7I < 0.05. The appropriate fre- 
quency resolution A/ is that corresponding to the 90-sec 
segment discrete Fourier transforms, so A/ = 0.011 Hz. 
As described later in Sec. IVII we calculate a statistical 
error, cto, associated with the stochastic background es- 
timate Yopt/Tscg,- Under the implicit assumption made 
in Eg. 15.231 that the detector noise is stationary, one can 
show that cry = Tcrj^. Finally, referring to Table HVl for 
an estimate of a^, and using the total H2-L1 observation 
time of 51 hours, we obtain ay ~ 2.8 x 10^ sec. Thus, 
the condition of Eq. E23 becomes: [Ti2{fL)Y/'^ < 1. 

Using this example estimate guide, specific lines 
are rejected by calculating the coherence function be- 
tween detector pairs for the full sets of analyzed SI data, 
and eliminating any frequency bins at which ri2(/i) > 
10~2. The coherence functions are calculated with a 
frequency resolution of 0.033 Hz, and approximately 
20,000 (35,000) averages for the LHO-LLO (LHO-LHO) 
pairs, corresponding to statistical uncertainty levels crp = 
1/A^avg of approximately 5 x 10~^ (3 x 10~^). The exclu- 
sion threshold thus corresponds to a cut on the coherence 
data of order 100 crp. 

For the H2-L1 pair, this procedure results in eliminat- 
ing the 250 Hz frequency bin, whose coherence level was 
about 0.02; the H2-L1 coherence function over the anal- 
ysis band is shown in Fig. [7| For H1-H2, the bins at 
168.25 Hz and 168.5 Hz were eliminated, where the co- 
herence was also about 0.02 (see ¥ig.^. The sources of 
these lines are unknown. For Hl-Ll, no additional fre- 
quencies were removed by the coherence threshold (see 
Fig.©. 
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FIG. 7: Coherence between the H2 and LI detector outputs 
during SI. The coherence is calculated with a frequency res- 
olution of 0.033 Hz and A'^avg ~ 20, 000 periodogram averages 
(50% overlap); Hann windows are used in the Fourier trans- 
forms. There are significant peaks at harmonics of 16 Hz 
(data acquisition buffer rate) and at 250 Hz (unknown origin). 
These frequencies are all excluded from the cross-correlation 
sum. The broadband coherence level corresponds to the ex- 
pected statistical uncertainty level of 1/A^avg ~ 5 x 10~^. 



It is worth noting that correlations at the n • 60 Hz 
lines are suppressed even without explicitly eliminating 
these frequency bins from the sum. This is because these 
frequencies have a high signal-to-noise ratio in the power 
spectrum estimates, and thus they have relatively small 
values in the optimal filter. The optimal filter thus tends 
to suppress spectral lines that show up in the power spec- 
tra. This effect is illustrated in Fig. |51 and is essentially 
the result of having four powers of Si(/) in the denomi- 
nator of the integrand of the cross-correlation, but only 
two powers in the numerator. Nonetheless, we chose to 
remove the n ■ 60 Hz bins from the cross-correlation sum 
for robustness, and as good practice for future analyses, 
where improvements in the electronics instrumentation 
may reduce the power line coupling such that the opti- 
mal filter suppression is insufficient. 

Such optimal filter suppression does not occur, how- 
ever, for the 16 Hz line and its harmonics, and the ad- 
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FIG. 8; Coherence between the HI and LI detector outputs 
during SI, calculated as described in the caption of Fig. |7| 
There are significant peaks at harmonics of 16 Hz (data ac- 
quisition buffer rate). These frequencies are excluded from 
the cross-correlation sum. The broadband coherence level 
corresponds to the expected statistical uncertainty level of 
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ditional 168.25, 168.5, 250 Hz lines; these lines typically 
do not appear in the power spectrum estimates, or do so 
only with a small signal-to-noise ratio. These lines must 
be explicitly eliminated from the cross-correlation sum. 
These discrete frequency bins are all zeroed out in the 
optimal filter, so that each excluded frequency removes 
0.25 Hz of bandwidth from the calculation. 



D. Results and interpretation 

The primary goal of our analysis is to set an upper 
limit on the strength of a stochastic gravitational wave 
background. The cross-correlation measurement is, in 
principle, sensitive to a combination of a stochastic grav- 
itational background and instrumental noise that is cor- 
related between two detectors. In order to place an upper 
limit on a gravitational wave background, we must have 
confidence that instrumental correlations are not playing 
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FIG. 9: Power spectral densities and optimal filter for the Hl- 
H2 detector pair, using the sensitivities shown in Fig. A 
scale factor has been applied to the optimal filter for display 
purposes. Note that spectral lines at 60 Hz and harmonics 
produce corresponding deep 'notches' in the optimal filter. 



a significant role. Gaining such confidence for the corre- 
lation of the two LHO interferometers may be difficult, 
in general, as they are both exposed to many of the same 
environmental disturbances. In fact, for the SI analy- 
sis, a strong (negative) correlation was observed between 
the two Hanford interferometers, thus preventing us from 
setting an upper limit on ^loo using the II1-II2 pair 
results. The correlated instrumental noise sources, rela- 
tively broadband compared to the excised narrow-band 
features described in the previous section, produced a 
significant II1-II2 cross-correlation (signal-to-noise ratio 
of —8.8); see Sec. I Villi for further discussion of the II1-II2 
instrumental correlations. 

On the other hand, for the widely separated (LHO- 
LLO) interferometer pairs, there are only a few paths 
through which instrumental correlations could arise. 
Narrow-band inter-site correlations are seen, as described 
in the previous section, but the described measures have 
been taken to exclude them from the analysis. Seismic 
and acoustic noise in the several to tens of Hz band 
have characteristic coherence lengths of tens of meters 
or less, compared to the 3000 km LHO-LLO separation, 
and pose little problem. Globally correlated magnetic 
field fluctuations have been identified in the past as the 
most likely candidate capable of producing broadband 
correlated noise in the widely separated detectors [33| . 
An order-of-magnitude analysis of this effect was made 
in Ref. concluding that correlated field fluctuations 
during magnetically noisy periods (such as during thun- 
derstorms) could produce a LHO-LLO correlated strain 
signal corresponding to a stochastic gravitational back- 
ground floh^QQ of order 10~*. These estimates evalu- 
ated the forces produced on the test masses by the corre- 
lated magnetic fields, via magnets that are bonded to the 



test masses to provide position and orientation control. -"^^ 
Direct tests made on the LIGO interferometers indicate 
that the magnetic field coupling to the strain signal was 
generally much higher during SI — up to 10^ times greater 
for a single interferometer — than these force coupling es- 
timates. Nonetheless, the correspondingly modified es- 
timate of the equivalent due to correlated magnetic 
fields is still 5 orders-of-magnitude below our present sen- 
sitivity. Indeed, Figs. [3 and show no evidence of any 
significant broadband instrumental correlations in the SI 
data. We thus set upper limits on ^Iq h^QQ for both the 
Hl-Ll and the H2-L1 pair results. 

Accounting for the combination of a gravitational wave 
background and instrumental correlations, we define an 
effective Q, fioff, for which our measurement Vopt/Tsog 
provides an estimate: 

f^off hloO = ^pt/Tscg = {^0 + i^inst) ^IQO • (5-24) 

Note that r^inst (associated with instrumental correla- 
tions) may be either positive or negative, while Qq for 
the gravitational wave background must be non-negative. 
We calculate a standard two-sided, frequentist 90% con- 
fidence interval on fioff /iioo follows: 

[ficff /i?oo - l-65CTn,tot , f^off/iioo + 1-65 CTn, tot] (5.25) 

where trn.tot is the total estimated error of the cross- 
correlation measurement, as explained in Sec. IVII In a 
frequentist interpretation, this means that if the experi- 
ment were repeated many times, generating many values 
of rioff /iioo s-'^d (Tr2,tot, then the true value of flcs H^qq is 
expected to lie within 90 percent of these intervals. We 
establish such a confidence interval for each detector pair. 

For the Hl-Ll and H2-L1 detector pairs, we are con- 
fident in assuming that systematic broad-band instru- 
mental cross-correlations are insignificant, so the mea- 
surement of flcff h-iQQ is used to establish an upper limit 
on a stochastic gravitational wave background. Specifi- 
cally, assuming Gaussian statistics with fixed rms devia- 
tion, CTo.tot, we set a standard 90% confidence level upper 
limit on H^qq . Since the actual value of fig must be 
non-negative, we set the upper limit to 1.28 ctji. tot if the 
measured value of fles /i^qq is less than zero.^"'^ Explicitly, 

^^0 ^100 < max{r2off hf^^ , 0} + 1.28CTj2,tot ■ (5.26) 

Table ^ summarizes the results obtained in applying 
the cross-correlation analysis to the LIGO SI data. The 



The actual limit on Qq /Ijqq that appears in Ref. llSl is 10~^ , 
since the authors assumed a magnetic dipole moment of the test 
mass magnets that is a factor of 10 higher than what is actually 
used. 

We are assuming here that a negative value of f!(,fl qo '^'^'^ 
random statistical fluctuations in the detector outputs and not 
to systematic instrumental correlations. 
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most constraining (i.e., smallest) upper limit on a gravi- 
tational wave stochastic background comes from the H2- 
Ll detector pair, giving Oq H^qq < 23. The significant 
H1-H2 instrumental correlation is discussed further in 
Sec. IVIIII The upper limits in Table ^ can be com- 
pared with the expectations given in Fig. ^ properly 
scaling the latter for the actual observation times. The 
H2-L1 expected upper limit for 50 hours of data would 
be flo /ifpQ < 14. The difference between this number 
and our result of 23 is due to the fact that, on average, 
the detector strain sensitivities were poorer than those 
shown in Fig. ^ 

In computing the Table ITTl numbers, some data selec- 
tion has been performed to remove times of higher than 
average detector noise. Specifically, the theoretical vari- 
ances of all 900-sec intervals are calculated, and the sum 
of the ay^^ is computed. We then select the set of largest 

CTy^ (i.e., the most sensitive intervals) that accumulate 
95^ of the sum of all the weighting factors, and include 
only these intervals in the Table UTI results. This selection 
includes 75 — 85% of the analyzed data, depending on the 
detector pair. We also excluded an additional ^ 10 hours 
of H2 data near the beginning of SI because of large data 
acquisition system timing errors during this period. The 
deficits between the observation times given in Table UTI 
and the total SI double-coincident times given in Sec. UTI 
are due to a combination of these and other selections, 
spelled out in Table ITTTI 

Shown in Fig. arc the weighting factors ay^^ (cf. 



Eq. 



over the duration of the SI run. The enter 



the expression for Y^pt (Eq. I5.13|l and give a quantita- 
tive measure of the sensitivity of a pair of detectors to a 
stochastic gravitational wave background during the I- 
th interval. Additionally, to gauge the accuracy of the 
weighting factors, we compared the theoretical standard 
deviations , to the measured standard deviations sy^ j 
(cf. Eg. I5.12|l . For each interferometer pair, all but one 
or two of the ctyu/syu ratios lie between 0.5 and 2, and 
show no systematic trend above or below unity. 

Finally, Fig. ^2 shows the distribution of the cross- 
correlation values with mean removed and normalized by 
the theoretical standard deviations — i.e., xjj = {Yjj — 
Yopt)/<7Yij- The values follow quite closely the expected 
Gaussian distributions. 
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FIG. 10: The weighting factors o-y^^ for each interferome- 
ter pair over the course of the SI run; each point represents 
900-sec of data. In each plot, a horizontal line indicates the 
weighting factor corresponding to the detector power spectra 
averaged over the whole run. 



can be quantitatively compared to the theoretically ex- 
pected signal arising from a stochastic background with 
^gw{f) = f^o = const by forming the statistic: 



E. Frequency- and time-domain characterization 

Because of the broadband nature of the interferome- 
ter strain data, it is possible to explore the frequency- 
domain character of the cross-correlations. In the anal- 
ysis pipeline, we keep track of the individual frequency 
bins that contribute to each Yjj, and form the weighted 
sum of frequency bins over the full processed data to pro- 
duce an aggregate cross-correlation spectrum, Yopt [I] , for 
each detector pair (cf. Eg. I5.16|l . These spectra, along 
with their cumulative sums, arc shown in Fig. 1121 l^pt[^] 



[Re{Yopt[e]) - Hy^^ 



(5.27) 



which is a quadratic function of CIq. The sum runs over 
the ~ 1000 frequency bins^^ contained in each spectra. 
The expected values /iy ^ [£] and theoretical variances 



To be exact, 1020 frequency bins were used for the Hl-Ll, H2-L1 
correlations and 1075 bins for H1-H2. 
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Interferometer fleg ftfoo 


f2eff hiQo/aQ^tot 


90% confidence 


90% confidence 


Xmin 


Frequency 


Observation 


pair 




interval on QcS ftioo 


upper limit 


(per dof) 


range 


time 


H1-H2 -8.3 


-8.8 


[-9.9 ±2.0, -6.8 ± 1.4] 




4.9 


40 - 300 Hz 


100.25 hr 


Hl-Ll 32 


1.8 


[2.1 ±.42, 61 ±12] 


fio hloo < 55 ± 11 


0.96 


40 - 314 Hz 


64 hr 


H2-L1 0.16 


0.0094 


[-30 ±6.0, 30 ±6.0] 


S^o hloo < 23 ± 4.6 


1.0 


40 - 314 Hz 


51.25 hr 



TABLE H: Measured 90% confidence intervals and upper limits for the three LIGO interferometer pairs, assuming fig„(/) = 
Qo = const in the specified frequency band. For all three pairs we compute a confidence interval according to Eq. 15.251 For 
the LHO-LLO pairs, we are confident in assuming the instrumental correlations are insignificant, and an upper limit on a 
stochastic gravitational background is computed according to Eq. 15.261 Our established upper limit comes from the H2-L1 
pair. The ± error bars given for the confidence intervals and upper limit values derive from a ±10% uncertainty in the 
calibration magnitude of each detector; see Sec. IVII and Table IIVI The Xmin per degree of freedom values are the result of a 
frequency-domain comparison between the measured and theoretically expected cross-correlations, described in Sec. IV El 



Selection criteria 


H1-H2 


Hl-Ll 


H2-L1 


A: All double- 


188 hr 


116 hr 


131 hr 


coincidence data 


46% 


28% 


32% 


B: A plus Tiock > 900-sec, 


134 hr 


75 hr 


81 hr 


& calibration monitored 


33% 


18% 


20% 


C: B plus valid GPS timing. 


119 hr 


75 hr 


66 hr 


& calibrations within bounds 


29% 


18% 


16% 


D: C plus quietest 


100 hr 


64 hr 


51 hr 


intervals 


25% 


16% 


13% 



Note that by using Eas. 15.71 15.81 one can show 



TABLE III: Summary of the selection criteria applied to the 
double-coincidence data for SI. A: portion of the 408- hr SI 
run having double-coincidence stretches greater than 600-sec; 
B: data portion satisfying criterion A, plus: data length is 
> 900-sec for the analysis pipeline, and the calibration mon- 
itoring was operational; C: data portion satisfying criterion 
B, plus: GPS timing is valid and calibration data are within 
bounds; D: data portion satisfying criterion C, plus: qui- 
etest data intervals that accumulate 95% of the sum of the 
weighting factors. This last data set was used in the analysis 
pipeline. 



CT~ [P\ are given by 
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vifhich are the expected value and theoretical variance of 
the weighted cross-correlation Fopt (c-f- Eg. I5.13|l . 

For each detector pair, we find that the minimum 
value occurs at the corresponding estimate fieff /iioo f^^' 
that pair, and the width of the = Xmin ^ 2. 71 points 
corresponds to the 90% confidence intervals given in Ta- 
ble |nj For the Hl-Ll and H2-L1 pairs, the minimum 
values are Xmin = (0.96, 1.0) per degree-of-freedom. This 
results from the low signal-to-noise ratio of the measure- 
ments: Oeff /ifoo/o'o^tot = (1.8, .0094). 

For the H1-H2 pair, Xmin — 4-9 per degree-of-frecdom. 
In this case the magnitude of the cross-correlation signal- 
to-noise ratio is relatively high, flgs /ifoo/'^'fJ-tot = —8.8, 
and the value of Xmin indicates the very low likelihood 
that the measurement is consistent with the stochastic 
background model. For ~I000 frequency bins (the num- 
ber of degrees-of-freedom of the fit), the probability of 
obtaining such a high value of Xmin i^ extremely small, 
indicating that the source of the observed H1-H2 correla- 
tion is not consistent with a stochastic gravitational wave 
background having f2gw(/) = Qq ~ const. 

It is also interesting to examine how the cross- 
correlation behaves as a function of the volume of data 
analyzed. Figure^lshows the weighted cross-correlation 
statistic values versus time, and the evolution of the esti- 
mate r^off ft-ioo Slid statistical error bar over the data 
run. Also plotted are the probabilities p(|rieff/iiool — 
I a; I) = 1 — erfd^l/v^ (To) of obtaining a value of r^cff ft^ioo 
greater than or equal to the observed value, assuming 
that these values are drawn from a zero-mean Gaussian 
random distribution, of width equal to the cumulative 
statistical error at each point in time. For the Hl-Ll 
and H2-L1 detector pairs, the probabilities are > 10% 
for the majority of the run. For H1-H2, the probability 



17 




100 150 200 

Frequency (Hz) 



250 300 



FIG. 12: Real part of the cross-correlation spectrum, 
for all I, J included in the TableEl ^/^-f (7^*^°^ Hz" ), for each detector pair. The grey 



FIG. 11: Normal probabilities and histograms of the values 

XIJ = (Ylj — Yopt)/(TYi_ 

results. In theory, these values should be drawn from a Gaus- 
sian distribution of zero mean and unit variance. The solid 
lines indicate the Gaussian that best fits the data; in the cu- 
mulative probability plots, curvature away from the straight 
lines is a sign of non-Gaussian statistics. 



line in each plot shows the cumulative sum of the spectrum 
from 40 Hz to ft, multiplied by 5f (which makes it dimen- 
sionless); the value of this curve at the right end is our es- 
timate 17cff /iioo- Note that the excursions in the cumulative 
sum for the Hl-Ll and H2-L1 correlations have magnitudes 
less than 1-2 error bars of the corresponding point estimates; 
simulations with uncorrelated data show the same qualitative 
behavior. 



drops below 10 after about 11 days, suggesting the 
presence of a non-zero instrumental correlation (see also 
Sec.lVIll. 
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FIG. 13: Cross-correlation statistics as a function of amount of data analyzed. Each column of plots shows the analysis results 
for a given detector pair, as indicated, over the duration of the data set. Top plots: Points correspond to the cross-correlation 
statistic values Yij appropriately normalized, MT^^ay^j Yi j / j ay'^^ , where M is the total number of analyzed intervals, 

so that the mean of all the values is the final point estimate f2eff /iioo- The scatter shows the variation in the point estimates 
from segment to segment. Middle plots: Evolution over time of the estimated value of fieff /ifoo- The black points give the 
estimates f2cff ^loo ^-^d the grey points give the ±1.65 an errors (90% confidence bounds), where an is defined by Eg. 16.11 The 
errors decrease with time, as expected from a T~^/^ dependence on observation time. Bottom plots: Assuming that the 
estimates shown in the middle plots are drawn from zero-mean Gaussian random variables with the error bars indicated, the 
probability of obtaining a value of \QcS /iiool ^ observed absolute value is given by: p{\QcH /iiool ^ l^^l) = 1 — erf(|j;|/\/2 an). 
This is plotted in the bottom plots. For the H1-H2 pair, the probability becomes < 10"^" after approximately 11 days. While 
the Hl-Ll pair shows a signal-to-noise ratio above unity, ftioo/S'fi,tot = 1.8, there is a 10% probability of obtaining an equal 
or larger value from random noise alone. 



F. Time shift analysis 

It is instructive to examine the behavior of the cross- 
correlation as a function of a relative time shift r intro- 
duced between the two data streams: 

Y{t) = / dh dt2Si{h-T)Q{ti^t2)s2{t2) 

J-T/2 J-T/2 

/oo 
df e^^^f^YU), (5.32) 
OO 

where 

/oo 
df stu^ nK{f)Q{ns2{n- (5.33) 
-00 



Thus, Y{t) is simply the inverse Fourier transform of the 
integrand, !"(/), of the cross-correlation statistic Y (c.f. 
Eq. I4.3|l . The discrete frequency version of this quan- 
tity, Vopti^] (c.f. Eg. I5.16|) . is shown in Fig. ^] for each 
detector pair. Figure IT^ shows the result of performing 
discrete inverse Fourier transforms on these spectra. For 
time shifts very small compared to the original FFT data 
length of 90-sec, this is equivalent to shifting the data and 
recalculating the point estimates. 

Also shown are expected time shift curves in the 
presence of a significant stochastic background with 
r2g„(/) = rig = const. These are obtained by taking 
the inverse discrete Fourier transforms of Eq. 15.281 they 
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have an oscillating behavior reminiscent of a sine func- 
tion. For the LHO-LLO pairs, these are computed for 
the upper limit levels given in Table while for Hl- 
H2, the expected curve is computed taking f2o ^loo equal 
to the instrumental correlation level of —8.3. For the 
two inter-site correlations (Hl-Ll and H2-L1), most of 
the points lie within the respective standard error lev- 
els: (Tfi^tot = 18 for Hl-Ll, and (7n,tot = 18 for H2-L1; 
for H1-H2, most of the points lie outside the error level, 
CTf2.tot = 0.95, indicating once again that the observed 
H1-H2 correlation is inconsistent with the presence of a 
stochastic background with rig„(/) = Qq = const. 



• H1-L1 data 
- - Expected for n = 55 




-80 

-0.04 -0.03 -0.02 -0.01 0.01 
Time shift (sec) 



0.02 0.03 



0.04 




-60 

-0.04 -0.03 -0.02 -0.01 0.01 0.02 0.03 0.04 
Time shift (sec) 

10 



^ 



• H1-H2data T 
- - Expected for £1 = -8.3 r f " 



-0°.04 -0.03 -0.02 -0.01 0.01 0.02 0.03 0.04 
Time shift (sec) 



VI. ERROR ESTIMATION 

We have identified three potentially significant types of 
error that contribute to the total error on our estimate 
of r^off h^QQ. The first is a theoretical statistical error, ^"^ 



(6.1) 



where the second equality follows from the definition of 
^cS ^100 ill terms of l^pt and 1/, treating the weighting 
factors as constants in the calculation of the theoret- 
ical variance of l^pt . We estimate this error by replacing 
the theoretical variance cryj(= (Ty^^/10) by its unbiased 
estimator Sy^^/10 (c.f. Eas. IsTtI 1^1^ .^^ Thus, 
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The last two sources of error are due to unresolved 
time variations in the interferometers' calibration, crn.cah 
and strain noise power spectra, crn,psd- As described 
earlier, detector power spectrum estimates are made on 
900-sec data intervals, and a single response function, de- 
rived from the central 60-seconds of calibration line data, 
is applied to each interval. Variations in both the re- 
sponse functions and the power spectra occur on shorter 
time scales, and we have estimated the systematic er- 
rors (an, cat and aa.psd) due to these variations as fol- 
lows. The cross-correlation analysis is performed again 
using a finer time resolution for calibration and power 
spectrum estimation, and the results are assumed to be 
representative of the effect of variations at other time 
scales. Specifically, each detector pair is re-analyzed with 
power spectrum estimates, and corresponding optimal fil- 
ters, computed for each 90-sec data segment (using the 
same frequency resolution as the original analysis, but 
approximately 1/10 the number of averages). Separately, 
each detector pair is also re-analyzed using the calibra- 
tion line amplitudes, and resulting response functions, 
corresponding to each 90-sec data segment. Each anal- 
ysis yields a new point estimate fioff/iioo! for each re- 
analysis, the difference between the new point estimate 
and the original point estimate is used as the estimate of 
the systematic errors trn.cai and an.psd- The total error 
is then formed as: 
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FIG. 14: Results of a time shift analysis for the three detector 
pairs. Plotted are the discrete inverse Fourier transforms of 
Yipt[Z]. Also shown are the expected time shift curves in the 
presence of a stochastic background with ngw(/) = Q,q = 
const; for the Hl-Ll and H2-L1 pairs, these are computed 
using the corresponding upper limit levels Slo ^iioo = 55 and 
^0 Ziioo = 23, respectively, while for H1-H2 the instrumental 
correlation level of —8.3 is used. 



Here we are treating Hoff loo ' ^opt , and 1/ as random variables 
and not as their values for a particular realization of the data. 
This is valid provided the individual cross-correlation measure- 
ments Yjj are statistically independent of one another. This 
assumption was tested by computing the auto-covariance func- 
tion of the Y[j data sequences; for each of the three Yjj sets, 
the result was a delta function at zero-lag, as expected for inde- 
pendent data samples. 
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These error estimates are shown in Table IIVI Also 
shown in the table are values for the fractional calibration 
uncertainty. The significant effect here is a frequency- 
independent uncertainty in the response function magni- 
tude; uncertainties in the phase response are negligible 
in the analysis band. 

We have also considered the effect of data acquisition 
system timing errors on the analysis. The behavior of the 
cross-correlation statistic when a time offset is introduced 
into the analysis was shown in Fig. El A finer resolution 
plot of the time-shift curve in the presence of a significant 
stochastic background indicates that a, t — ±400 /isec 
offset between the LHO and LLO interferometers corre- 
sponds to a 10% reduction in the estimate of our upper 
limit. The growth of this error is roughly quadratic in 
T. Throughout the SI run. the time-stamping of each in- 
terferometer's data was monitored, relative to GPS time. 
The relative timing error between H2 and LI was approx- 
imately 40 fiscc for roughly half the analyzed data set, 
320 /isec for 32% of the data set, and 600 /isec for 16% 
of the set. The combined effect of these timing offsets 
is an effective reduction in the point estimate, ftcs hf^Q, 
of 3.5%, a negligible effect. The Hl-Ll relative timing 
errors were even smaller, being less than 30 ^sec during 
the whole data set. 



VII. VALIDATION: SIGNAL INJECTIONS 

The analysis pipeline was validated by demonstrating 
the ability to detect coherent excitation of the interferom- 
eter pairs produced by simulated signals corresponding 
to a stationary, isotropic stochastic gravitational wave 
background. 

A software package was developed to generate a pseu- 
dorandom time series representing this excitation. In 
this manner, pairs of coherent data trains of simulated 
stochastic signals could be generated. The amplitude of 
the simulated stochastic background signal was adjusted 
by an overall scale factor, and the behavior of the detec- 
tion algorithm could be studied as a function of signal- 
to-noise ratio. Simulated data were either injected into 
the interferometer servo control system in order to di- 
rectly stimulate the motion of the interferometer mirrors 
(hardware injections) or the calculated waveforms could 
be added in software to the interferometer data as part 
of the analysis pipeline (software injections). The former 
approach was used to inject a few simulated stochastic 
background signals of different amplitudes during inter- 
ferometer calibrations at the beginning and end of the SI 
run. The latter approach was used after the SI run dur- 
ing the data analysis phase. Table lists the different 
injections that were used to validate our procedure. 



A. Hardware injected signals 

Hardware injections required that the simulated data 
trains be first convolved with the appropriate instrument 
response functions. These pre-processed data trains were 
then injected digitally into the respective interferometer 
servo control systems. 

Simulated stochastic background signals with 
figw(/) = flo — const were injected simultaneously 
into the H2-L1 pair for two 1024 second (17.07 min) 
periods shortly after the SI run was completed. Refer- 
ring to Table the two injections had different signal 
strengths, corresponding to signal-to-noise ratios of ^ 20 
and ~ 10, respectively. The stronger injection produced 
a noticeable increase in the H2 power spectrum in the 
band from 40 to 600 Hz. 

In principle, the stochastic gravitational wave back- 
ground estimate can be derived from a single (point) 
measurement of the cross-correlation between pairs of in- 
terferometers. However, in order to verify that the sim- 
ulated signals being detected were consistent with the 
process being injected, time shift analyses of the data 
streams were performed for a number of different offsets, 
T. This technique can potentially identify instrumental 
and environmental correlations that are not astrophysi- 
cal. 

For each injection, the results of the time-shifted anal- 
ysis were compared with the expected time shift curves. 
Allowing for possible (unknown) time shifts associated 
with the stimulation and data acquisition processes, a 
two parameter regression analysis was performed on 
the time shift data to determine: (i) the time offset (if 
any existed), (ii) the amplitude of the signal, and (iii) the 
uncertainties in the estimation of these parameters. 

Results of this analysis for the hardware injection with 
signal strength JIq /if qq = 3906 are shown in Fig. E| The 
agreement between injected simulated signal and the de- 
tected signal after end-to-end analysis with our pipeline 
gave us confidence that the full data analysis pipeline was 
working as expected. 



B. Software injected signals 

The same simulated signals can be written to file, and 
then added to the interferometer strain channels. These 
software simulation signals were added in after the strain 
data were decimated to 1024 Hz, as shown in Fig. O 
The flexibility of software injection allowed a wide range 
of values for ^Iq H^qq to be studied. Refer to Table |Vl 
for details. This allowed us to follow the performance 
of the pipeline to smaller signal-to-noise ratios, until the 
signal could no longer be distinguished from the noise. 
The behavior of the deduced signal versus injected signal 
at a large range of signal-to-noise ratios is presented in 
Fig. HI 
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Pair 


an 






$S2,tot 


Calibration 
uncertainty 


H1-H2 


0.93 


0.078 


0.16 


0.95 


±20% 


Hl-Ll 


18 


0.23 


0.29 


18 


±20% 


H2-L1 


15 


9.3 


1.2 


18 


±20% 



TABLE IV: Sources of error in the estimate QcS ^loo = Kjpt/Tseg: o^n is the statistical error; Sn.psd is the error due to unresolved 
time variations of the equivalent strain noise in the detectors; and an,cai is the error due to unresolved calibration variations. 
The calibration uncertainty for each detector pair results from adding linearly a ±10% uncertainty for each detector, to allow 
for a worst case combination of systematic errors. 
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Magnitude of Detected Signal 


Pair 


Software (SW) 
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SNR 


(90% CL, Qo ftMo±l-65Sn) 


H2-L1 


HW 


3906 
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3744 ± 663 


H2-L1 


HW 


24414 


17 


25365 ± 2341 


H2-L1 


SW 


16 






H2-L1 


SW 


100 






H2-L1 


SW 


625 


3 


891 ± 338 


H2-L1 


SW 


3906 


13 


4361 ± 514 


H2-L1 


SW 


24414 


50 


25124 ± 817 



TABLE V: Summary of injected signals used to validate the analysis pipeline. Both hardware injections during the SI run and 
post-Sl software injections were used. Injections were introduced into short data segments (refer to text). The signal-to-noise 
ratios shown correspond to integration times that are much shorter than the full SI data set, and thus the lower signal-to-noise 
ratio injections were not detectable. The software and hardware injections have different signal-to-noise ratios for the same 
^0 ^100 values due to the variation in the interferometer noise power spectral densities at the different epochs when the signals 
were injected. 



VIII. THE H1-H2 CORRELATION 

The significant instrumental correlation seen between 
the two LHO interferometers (HI and H2) prevents us 
from establishing an upper limit on the gravitational 
wave stochastic background using what is, potentially, 
the most sensitive detector pair. It is thus worth exam- 
ining this correlation further to understand its charac- 
ter. We tested the analysis pipeline for contamination 
from correlated spectral leakage by re-analyzing the Hl- 
H2 data using a Hann window on the 90-sec data seg- 
ments instead of the Tukey window. The result of this 
analysis, when scaled for the effective reduction in obser- 
vation time, was — within statistical error — the same as 
the original Tukey-windowed analysis, discounting this 
hypothesis. 

Some likely sources of instrumental correlations are: 
acoustic noise coupling to both detectors through the 
readout hardware (those components not located in- 
side the vacuum system) ; common low- frequency seismic 
noise that bilinearly mixes with the 60 Hz and harmonic 
components, to spill into the analysis band. Figure 1171 
shows the coherence function (Eg. 15.2^ between HI and 
H2, calculated over approximately 150 hours of coinci- 
dent data. It shows signs of both of these types of sources. 

Both of these noise sources are addressable at the in- 



strument level. Improved electronics equipment being 
implemented on all detectors should substantially reduce 
the n ■ 60 Hz lines, and consequently the bilinearly mixed 
sideband components as well. Better acoustic isolation 
and control of acoustic sources is also being planned 
to reduce this noise source. It is also conceivable that 
signal processing techniques, such as those described in 
Ref. HI Hi, could be used to remove correlated noise, 
induced by measurable environmental disturbances, from 
the data. 



IX. CONCLUSIONS AND FUTURE PLANS 

In summary, we have analyzed the first LIGO science 
data to set an improved, direct observational upper limit 
on a stochastic background of gravitational waves. Our 
90% confidence upper limit on a stochastic background, 
having a constant energy density per logarithmic fre- 
quency interval, is flo /if qq < 23 in the frequency band 
40 — 314 Hz. This is a roughly 10* times improvement 
over the previous, broadband interferometric detector 
measurement. 

We described in detail the data analysis pipeline, and 
tests of the pipeline using hardware and software injected 
signals. We intend to use this pipeline on future LIGO 
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FIG. 15: Hardware injection time shift analysis for the H2- 
Ll interferometer pair, with signal strength Q,o /iioo = 3906. 
Panel (a) : Time shift dependence of the cross-correlation (re- 
fer to Eg. 15.321 . The data are shown with ±lan error bars 
estimated from the measured quantities (Eq. 16.11 for each 
time offset. The dashed curve is the expected dependence, 
scaled and offset in time to provide a best fit. Panel (b): 
Contour plot of x^i^o,''') near the best fit. The minimum 
value is Xmin ~ 1-8 (for 2 degrees of freedom), and occurs 
at the coordinates of the black rectangle: {Slo ftioO) ~ 
{3744, —270 usee}. The cross (-I-) corresponds to the in- 
jected signal, whose estimated strength has 90% confidence 
bounds of: 3345 < Qo ^loo < 4142. The best fit time oflFset 
of —270 fisec is within the observed relative data acquisition 
timing errors between H2 and LI during SI. 
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FIG. 16: H2-L1 point estimates and error bars obtained from 
the SI data analysis for both hardware and software injec- 
tions. Measured versus injected SNRs are shown for a num- 
ber of simulations. The ordinate of each point is the result of 
a analysis like the one shown in Fig. 1151 The x? fit a-lso 
provides an estimate an of the measurement noise. The esti- 
mate is used to normalize the measured and injected values 
of f2() . 



an upper limit on JIq /iiqo • Eventually, with both 4-km 
interferometers (HI and LI) operating at the design sen- 
sitivity level shown in Fig. ^ we expect to be able to set 
an upper limit using 1-year of data from this detector 
pair at a level f^o /i?oo < 1 x lO"*" in the 40 - 314 Hz 
band. This would improve on the limit from big-bang 
nucleosynthesis (see Table IJ). The two interferometers 
at LHO (HI and H2) could potentially provide a lower 
upper limit, but given our present level of correlated in- 
strumental noise (jilinst /iiool ^ 10) j first need to re- 
duce the correlated noise in each detector by a factor of 
^ 10"*. We also anticipate cross-correlating LI with the 
ALLEGRO resonant bar detector (located nearby LLO 
in Baton Rouge, LA) for a higher frequency search. With 
this pair performing at design sensitivity, an upper limit 
of order h\^Q < 0.01 could be set around 900 Hz, using 
1 year of coincident data. 



science data to set upper limits on fio h^QQ at levels which 
are orders of magnitude below unity. Two possible addi- 
tions to the treatment presented here are being consid- 
ered for future analyses: a method for combining upper 
limits from Hl-Ll and H2-L1 that takes into account the 
potential H1-H2 instrumental correlations; a Bayesian 
statistical analysis for converting the point estimate into 
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puts during SI. The coherence is calculated with a frequency 
resolution of 0.033 Hz, and with approximately 35,000 peri- 
odogram averages; 50% overlap Hann windows were used in 
the Fourier transforms. In addition to the near-unity coher- 
ence at 60 Hz and harmonics, there is broadened coherence 
at some of these lines due to bilinear mixing of low-frequency 
seismic noise. There are several few-Hz wide regions of signif- 
icant coherence (around 168 Hz, e.g.), and the broad region 
of significant coherence between 220 Hz and 240 Hz, that 
are likely due to acoustic noise coupling. Also discernible 
are small peaks at many of the integer frequencies between 
245-310 Hz, likely due to coupling from the GPS 1 pulse-per- 
second timing signals. 

I 

APPENDIX A: LIST OF SYMBOLS 



Symbol 


Description 


E q. No. 


ng„(/) 


Energy density in gravitational waves per logarithmic frequency interval 






in units of the closure energy density pc 




Pc Pgw 


Critical energy density needed to close the universe, and total energy den- 


iniAi 




sity in gravitational waves 




/lloo, Ho, Hioo 


/iioo is the Hubble constant, Hq, in units of i^ioo = 100 km/sec/Mpc 




habit), h{t) 


Perturbations of the space-time metric, and the corresponding gravita- 






tional wave strain in a detector 




xo, x\ 


Position vector of an interferometer vertex, and unit vectors pointing in 


IlllAI 




the directions of the arms of an interferometer 




'S'gw(/) 


Power spectrum of the gravitational wave strain h(t) 


EH 



fifaiftttiillJi|;<ihJ>JhiiM|Mihli4iMYJi.i^u.AJi | lAte 



T r 

40-71 Hz 



aA4itw<iik|fc/1u»ifcii|iii^<<MLl,nki|<HjitoJ^iJ^^^ 
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- 190-221 Hz - 
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FIG. 17: Coherence between the HI and H2 detector out- 
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rib 



^min; -t-max 



Si{t), Si{f) Equivalent strain output of the i-th detector 14.11 

hiit), hi{f) Gravitational wave strain in the i-th detector 14.11 

Tiiit), ni{f) Equivalent strain noise in the i-th detector 14.11 

ri{t), ri[k] Raw (i.e., uncalibrated) output of the i-th detector for continuous and 15.11 
discrete time 

Ri{f), Rii[£] Response function for the i-th detector, and the discrete frequency re- IV Al 

spouse function for interval / 

tfe , fi Discrete time and discrete frequency values I V Al 

5t, 5f Sampling period of the time-series data (1/1024 sec after down-sampling), IV Al 

and bin spacing (0.25 Hz) of the discrete power spectra, optimal filter, . . . 

Af General frequency resolution of discrete Fourier transformed data IV Gl 

Number of discrete-time data points in one segment of data IV Al 

riij[k], giij[k], gijK[q] Raw detector output for the J-th segment in interval I evaluated at dis- 15.21 
Crete time tk , and the corresponding windowed and zero-padded time series 
and discrete Fourier transform 

Qij[P\ Cross-spectrum of the windowed and zero-padded raw time-series, binned 15.41 
to match the frequency resolution of the optimal filter Qi [(\ 

Window function for the i-th detector IV Al 

Number of frequency values binned together to match the frequency res- IV Al 
olution of the optimal filter Qi[P\ 

Indices corresponding to the maximum and minimum frequencies used in IV Al 
the calculation of the cross-correlation Yij 

St if) Finite-time approximation to the Dirac delta function S{f) 14.41 

T, Tscg, Tint General observation time, and durations of an individual data segment IV Al 

and interval (90 sec, 900 sec) 

Y General cross correlation of two detectors 14.21 

Q{t), Q{f) Optimal filter for the cross-correlation Y 14.81 

fj,Y, cry, Py Theoretical mean, variance, and signal-to-noise ratio of the cross- 14.514.61 
correlation Y 

(Py) Expected value of the signal-to-noise ratio of the cross-correlation Y 14.111 

Yij, Yi Gross-correlation for the J-th segment in interval /, and average of the 15.51 15. ifl 

Yij 

Yopt Weighted average of the Yj 15.131 

xij Gross-correlation values Yjj with mean removed and normalized by the IV Dl 
theoretical variances 

Yij[£], Yi[e], Yopt[e] Summands of Yij, Yj, Fopt ^MEMEM 

Qi[P\ Optimal filter for the cross-correlation Yjj 15.61 

7(/)i 7M Overlap reduction function evaluated at frequency /, and discrete fre- II VI IV Al 
quency fi 

7rms, Abw Root-mean-square value of the overlap reduction function over the corre- IIVI 

sponding frequency bandwidth Abw 

Pi{f), Pii[l] Power spectrum of the strain noise in the i-th detector, and the discrete IIVI IV Al 

frequency strain noise power spectrum estimate for interval / 

A/", A// Normalization factors for the optimal filter Q, Qi 14.91 15.71 

wTw2, 'wfw2 Overall multiplicative factors introduced by windowing 15.9115.1(71 

CTy , Sy Theoretical and estimated variance of the cross-correlation Yjj 15.81 15. 11?1 

ffYj 5 '^Yapt Theoretical variance of Yj and Vopt IVll 

/ir, [£] , cr|, [£] Theoretical mean and variance of Yopt [£] 15.281 15.291 

^opt Jopt ^^^^^ 

Pi2(/) Gross-power spectrum of the strain noise between two detectors IV Gl 

ri2(/) Goherence function between two detectors 15.221 

Ni^vg, o"r Number of averages used in the measurement of the coherence, and the IV Gl 

corresponding statistical uncertainity in the measurement 

Jlo, ^0 Actual and estimated values of an (assumed) constant value of r2gw(/) due IIII Al \W~K\ 

to gravitational waves 

fiinst, f^inst Actual and estimated values of the instrumental contribution to the mea- IV Dl 
sured cross-correlation 

QcH, f^cfi Actual and estimated values of an effective SI due to instrumental and 15.241 
gravitational wave effects 

aft, a'fi Actual and estimated variances of Q,cR due to statistical variations in Yopt 16.11 16.21 

ctq (,^1, S'q ^.j^i Actual and estimated variances of ficff due to variations in the instrument IVll 
calibration 
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2 ^2 


Actual and estimated variances oi i IcS due to variations m the noise power 


IV 11 


spectra 




^2 

""n.tot 


Estimated variance ot iZcfi due to combined statistical, calibration, and 


16.31 




puwci bpcCLia vaiiaLiulib 




x'(no) 


Chi-squared statistic to compare Fopt [£] to its expected value for a stochas- 






tic background with f2gw(/) = fio 




Xmin 


Minimum chi-squared value per degree of freedon 


EH 


^-W, Yif) 


Cross correlation statistic as a function of time-shift r, and its Fourier 


15.32115.331 



transform 



TABLE VI: A list of symbols that appear in the paper, along with their descriptions and equation numbers (if 
applicable) or sections in which they were defined. 



[1] B. Barish and R. Weiss, "LIGO and the Detection 
of Gravitational Waves," Phys. Today 52, 44 (1999). 
jhttp : / / www . li go . caltech. edu/ ^ 

[2] B. Willke et al., "The GEO 600 gravitational wave 
detector," Class. Quant. Grav. 19, 1377 (2002). 
|http://ww w . geo600 . iini-hannover Tde/| 

[3] B. Caron et al., "The VIRGO Interferometer for Grav- 
itational Wave Detection," Nucl. Phys. B-Proc. Suppl. 
54, 167 (1997). http://www.virgo.infn.it/ 

[4] K. Tsubono, "300-m laser interferometer gravitational 
wave detector (TAMA300) in Japan," in 1st Edoardo 
Amaldi Conf. on Gravitational Wave Experiments, edited 
by E. Coccia, G. Pizzella, and F. Ronga (Singapore, 
World Scientific, 1995), p. 112. 

[5] R. Savage, P. King, and S. Seel, "A highly stabilized 
10-watt Nd:YAG laser for the laser interferometer grav- 
itational wave observatory (LIGO)," Laser Phys. 8, 679 
(1998). 

[6] J. Giaime, P. Saha, D. Shoemaker, and L. Sievers, "A 
passive vibration isolation stack for LIGO: design, mod- 
eling, and testing," Rev. Sci. lustrum. 67, 208 (1996). 

[7] J.A. Giaime, E.J. Daw, M. Weitz, R. Adhikari, 
P. Fritschel, R. Abbott, R. Bork, and J. Heefner, "Feed- 
forward reduction of the microseism disturbance in a 
long-base-line interferometric gravitational-wave detec- 
tor," Rev. Sci. lustrum. 74, 218 (2003). 

[8] P. Fritschel, R. Bork, G. Gonzalez, N. Mavalvala, 
D. Ouimette, H. Rong, D. Sigg, and M. Zucker, "Read- 
out and control of a power-recycled interferometric 
gravitational-wave antenna," Appl. Opt. 40, 4988 (2001). 

[9] P. Fritschel, N. Mavalvala, D. Shoemaker, D. Sigg, 
M. Zucker, and G. Gonzalez, "Alignment of an inter- 
ferometric gravitational wave detector," Appl. Opt. 37, 
6734 (1998). 

[10] R. Adhikari, G. Gonzalez, M. Landry, and B. O'Reilly, 
"Calibration of the LIGO detectors for the First Scien- 
tific Run," in Proceedings of the 7th Gravitational Wave 
Data Analysis Workshop, Kizu-cho, Kyoto, Japan, 17-19 
December 2002 (unpublished). 

[11] The LIGO Scientific Collaboration, "Detector Descrip- 
tion and Performance for the First Coincidence Obser- 



vations between LIGO and GEO," gr-qc/0308043. Ac- 
cepted for publication in Nuclear Instruments and Meth- 
ods A. 

[12] B. Allen, "The stochastic gravity-wave background: 
sources and detection," in Proceedings of the Les Houches 
School on Astrophysical Sources of Gravitational Waves, 
Les Houches, 1995, edited by J. A. Marck and J. P. La- 
sota (Cambridge, 1996), p. 373. 

[13] M. Maggiore, "Gravitational wave experiments and early 
universe cosmology," Phys. Rep. 331, 283 (2000). 

[14] W. Freedman et al., "Final results from the Hubble Space 
Telescope Key project to measure the hubble constant," 
Astrophys. J. 553, 47 (2001). 

[15] B. K. Gibson and C. B. Brook, astro-ph/0011567; ASP 
Conference Series, 666 (2001). 

[16] C. L. Bennett et al., "First Year Wilkinson Microwave 
Anisotropy Probe (WMAP) Observations: Preliminary 
Maps and Basic Results," astro-ph/0302207. 

[17] C.W. Misner, K.S. Thorne, and J.A. Wheeler, Gravita- 
tion, (Freeman, San Francisco, 1973). 

[18] B. Allen and J.D. Romano, "Detecting a stochastic 
background of gravitational radiation: Signal processing 
strategies and sensitivities," Phys. Rev. D 59, 102001 
(1999). 

[19] R. R. Caldwell, R. A. Battye, E. P. S. SheUard, "Relic 
Gravitational Waves from Cosmic Strings: Updated Con- 
straints and Opportunities for Detection," Phys. Rev. D 
54, 7146 (1996). 

[20] D. M. Coward, R. R. Burman, and D. G. Blair, 
"The stochastic background of gravitational waves 
from neutron star formation at cosmological distances," 
Mon. Not. R. Astron. Soc. 324, 1015 (2001). 

[21] V. Ferrari, S. Matarrese, and R. Schneider, "Gravita- 
tional wave background from a cosmological population 
of core-collapse supernovae ," Mon. Not. R. Astron. Soc. 
303, 247 (1999). 

[22] B. Allen and A. C. Ottewill, "Detection of anisotropics 
in the gravitational-wave stochastic background," 
Phys. Rev. D 56, 545 (1997). 

[23] S. Drasco and E. E. Flanagan, "Detection methods 



26 



for non-Gaussian gravitational wave stochastic back- 
grounds," Phys. Rev. D 67, 082003 (2003). 

[24] M. S. Turner, "Detectability of infiation-produced grav- 
itational waves," Phys. Rev. D 55, R435 (1997). 

[25] M. P. McHugh, G. Zalamansky, F. Vernotte, E. Lantz, 
"Pulsar timing and upper limits on a gravitational wave 
background: A Bayesian approach," Phys. Rev. D 54, 
5993 (1996). 

[26] E. W. Kolb and M. Turner, The Early Universe, (Fron- 
tiers in Physics, Addison Wesley, 1990) . 

[27] J. Hough, J. R. Pugh, R. Bland, and R. W. Drever, 
"Search for continuous gravitational radiations," Nature 
254, 498 (1975). 

[28] K. Compton, D. Nicholson, and B. F. Schutz, in Proceed- 
ings of the Seventh Marcel Grossman Meeting on General 
Relativity, (World Scientific, 1994), p. 1078. 

[29] P. Astone et al., "Upper limit for a gravitational- wave 
stochastic background with the Explorer and Nautilus 
resonant detectors," Phys. Lett. B 385, 421 (1996). 

[30] P. Astone et al., "Upper limit at 1.8 kHz for a 
gravitational- wave stochastic background with the it AL- 
TAIR resonant-mass detector," Astronomy and Astro- 
physics 343, 19 (1999). 

[31] P. Astone et al., "Cross-correlation measurement of 
stochastic gravitational waves with two resonant grav- 
itational wave detectors," Astronomy and Astrophysics 
351, 811 (1999). 

[32] P. Astone, V. Ferrari, M. Maggiore, and J. D. Ro- 
mano, "Stochastic background of gravitational waves," 
Int. J. Mod. Phys. D 9, 361 (2000). 

[33] P.F. Michelson, "On detecting stochastic background 
gravitational radiation with terrestrial detectors," 
Mon. Not. R. Astron. Soc. 227, 933 (1987). 



[34] N. Christensen, "Measuring the stochastic gravitational 
radiation background with laser interferometric anten- 
nas," Phys. Rev. D 46, 5250 (1992). 

[35] E.E. Flanagan, "Sensitivity of the laser interferometer 
gravitational wave observatory (LIGO) to a stochastic 
background, and its dependence on the detector orienta- 
tions," Phys. Rev. D 48, 389 (1993). 

[36] The large-scale computing was carried out on Linux 
clusters at Caltech and the University of Wis- 
consin - Milwaukee (The UWM Medusa cluster, 
http: //www. Isc-group. phys .uwm. edu/beowulf /medusa), 
partially supported by the GriPhyN 
(http : / /www . griphyn . org) and iVDGL 

(http : / /www . ivdgl . org) projects. 

[37] W. H. Press, B. P. Flannery, S. A. Teukolsky, and 
W. T. Vetterling, in Numerical Recipes in C (Cambridge 
University Press, Cambridge, 1992), Second edition. 

[38] W. R. Johnston, J. D. Romano and J. T. Whelan, in 
preparation, 2003. 

[39] J.G. Proakis, D.G. Manolakis, in Digital Signal Process- 
ing: Principles, Algorithms, and Applications, (Prentice 
HaU, Upper Saddle River, N.J., 1996), Third edition, 
p. 624. 

[40] J. S. Bendat and A. G. Piersol, in Random Data; Anal- 
ysis and Measurement Procedures, (Wiley-Interscience, 
New York, N. Y., 2000), Third edition, p. 147. 

[41] B. Allen, W. Hua, and A. Ottewill, "Automatic cross- 
talk removal from multi-channel data," gr-qc/9909083. 

[42] A.C. Searle, S.M. Scott, and D.E. McClelland, "Spectral 
line removal in the LIGO data analysis system (LDAS)," 
Glass. Quant. Grav. (to appear), gr-qc/0308024. 



